Optimal network topologies for information transmission in active networks 
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This work clarifies the relation between network circuit (topology) and behavior (information 
transmission and synchronization) in active networks, e.g. neural networks. As an application, we 
show how to determine a network topology that is optimal for information transmission. By optimal, 
we mean that the network is able to transmit a large amount of information, it possesses a large 
number of communication channels, and it is robust under large variations of the network coupling 
configuration. This theoretical approach is general and does not depend on the particular dynamic 
of the elements forming the network, since the network topology can be determined by finding a 
Laplacian matrix (the matrix that describes the connections and the coupling strengths among the 
elements) whose eigenvalues satisfy some special conditions. To illustrate our ideas and theoretical 
approaches, we use neural networks of electrically connected chaotic Hindmarsh-Rose neurons. 
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The relation between neural circuits and behavior is 
a fundamental matter in neuroscience. In this work, 
we present a theoretical approach that has the poten- 
tial to unravel such a relationship in terms of network 
topology, information, and synchronization, in active net- 
works, networks formed by elements that are dynamical 
systems (such as neurons, chaotic or periodic oscillators). 
As a direct application of our approaches, we show how 
one can construct optimal neural networks that not only 
transmit large amounts of information from one element 
to another in the network, but also are robust under alter- 
ations in the coupling configuration. We also show that 
the relation between synchronization and information is 
rather subtle. Neural networks whose configurations al- 
low the transmission of large amounts of information 
might have at least two unstable modes of oscillation that 
are out of synchrony, while all the others are synchronous. 
Depending on the kind of measurement being done, one 
can arrive at contradicting statements concerning the re- 
lation between information and synchronization. We il- 
lustrate our theoretical approaches by using neural net- 
works of electrically connected chaotic Hindmarsh-Rose 
neurons [l|. These results have a tremendous impact in 
the understanding of information transmission in brain- 
like networks, as well as, in the mammalian brain. They 
also shed light on a better understanding of the neural 
code, the rules under which neurons encode and transmit 
information about external stimuli. 



II. BLURB 

This work shows how to relate in an active network 
the rate of information that can be transmitted from one 
point to another, regarded as mutual information rate 
(MIR), the synchronization level among elements, and 



III. INTRODUCTION 

Given an arbitrary time dependent stimulus that ex- 
ternally excites an active network formed by systems that 
have some intrinsic dynamics (e.g. neurons and oscilla- 
tors), how much information from such stimulus can be 
realized by measuring the time evolution of one of the 
elements of the network ? Determining how and how 
much information flows along anatomical brain paths is 
an important requirement for the understanding of how 
animals perceive their environment, learn and behave 

aii. 

The works of Refs. 0, i, i, i, i, 0| propose ways to 
quantify how and how much information from a stimulus 
is transmitted in neural networks. In particular, Ref. [5| 
demonstrated that 50% of the information about light 
displacements might be lost after being processed by the 
HI neuron, sensitive to image motion around a vertical 
axis, a neuron localized in a small neural network of the 
Chrysomya magacephala fly, the lobula plate. Does that 
mean that the HI neuron has an information capacity 
lower than the information contained in the light stim- 
ulus ? Or does that mean that information is lost due 
to the presence of internal noise ? These questions and 
others, which are still awaiting answers, concern the rules 
under which information is coded and then transmitted 
by neurons and it is a major topic of research in neuro- 
science referred to as the neural code [3, |j| . 

Even though the approaches of Ref. 3, i, i, H, i, 0| 
have brought considerable understanding on how and 
how much information from a stimulus is transmitted in 
a neural network, the relation between network circuits 
(topology) and information transmission in a neural as 
well as an active network is still awaiting a more quan- 
titative description [8"]. And that is the main thrust of 



the present manuscript, namely, to present a quantitative 
way to relate network topology with information in active 
networks. Since information might not always be easy to 
be measured or quantified in experiments, we endevour 
to clarify the relation between information and synchro- 
nization, a phenomenon! which is often not only possible 
to observe but also relatively easy to characterize. 

We initially proceed along the same line as in Refs. 
[3, Il3| J and study the information transfer in autonomous 
systems. However, instead of treating the informa- 
tion transfer between dynamical systems components, we 
treat the transfer of information per unit time exchanged 
between two elements in an autonomous chaotic active 
network. Thus, we neglect the complex relation between 
external stimulus and the network and show how to cal- 
culate an upper bound value for the mutual information 
rate (MIR) exchanged between two elements (a commu- 
nication channel) in an autonomous network. Ultimately, 
we discuss how to extend this formula to non-chaotic net- 
works suffering the influence of a time-dependent stimu- 
lus. 

Most of this work is directed to ensure the plausibil- 
ity and validity of the proposed formula for the upper 
bound of MIR (Sec. |V| and also to study its applications 
in order to clarify the relation among network topology, 
information, and synchronization. We do not rely only 
on results provided by this formula, but we also calcu- 
late the MIR by the methods in Refs. [ll|, d and by 
symbolic encoding the trajectory of the elements forming 
the network and then measuring the mutual information 
provided by this discrete sequence of symbols (method 
described in Sec. IXIII[) . 

To illustrate the power of the proposed formula, we ap- 
plied it to study the exchange of information in networks 
of coupled chaotic maps (Sec. IXVip and in Hindmarsh- 
Rose neural networks bidirectionally electrically coupled 
fSec lVI|) . The analyses are carried out using quantities 
that we believe to be relevant to the treatment of in- 
formation transmission in active networks: a communi- 
cation channel, the channel capacity, and the network 
capacity (see definitions in Sec. IXVp . 

A communication channel represents a pathway 
through which information is exchanged. In this work, 
a communication channel is considered to be formed by 
a pair of elements. One element represents a transmiter 
and the other a receiver, where the information about 
the transmiter can be measured. 

The channel capacity is defined in terms of the pro- 
posed upper bound for the MIR. It measures the local 
maximal rate of information that two elements in a given 
network are able to exchange, a point-to-point measure 
of information exchange. As we shall see, there are two 
network configurations for which the value of the upper 
bound can be considered to be maximal with respect to 
the coupling strength. 

The network capacity is the maximum of the KS- 
entropy, for many possible network configurations with 
a given number of elements. It gives the amount of in- 



dependent information that can be simultaneously trans- 
mitted within the whole network, and naturally bounds 
the value of the MIR in the channels, which concerns only 
the transmission of information between two elements. 

While the channel capacity is bounded and does not 
depend on the number of elements forming the network, 
the network capacity depends on the number of elements 
forming the network. 

As a direct application of the formula for the upper 
bound value of the MIR, we show that an active network 
can operate with a large amount of MIR and KS-entropy 
and at the same time it is robustly resistant to alterations 
in the coupling strengths, if the eigenvalues of the Lapla- 
cian matrix satisfy some specified conditions (Sec. IVIip . 
The Laplacian matrix describes the connections among 
the elements of the network. 

The conditions on the eigenvalues depend on whether 
the network is constructed in order to possess communi- 
cation channels that are either self-excitable or non-self- 
excitable (see definition in Sec. IXIVp . Active networks 
that possess non-self-excitable channels (formed by oscil- 
lators as the Rossler, or the Chua's circuit) have channels 
that achieve their capacity whenever their elements are in 
complete synchrony. Therefore, if a large amount of in- 
formation is desired to be transmitted point-to-point in 
a non-self-excitable network, easily synchronizable net- 
works are required. On the other hand, networks that 
possess self-excitable channels (as the ones formed by 
neurons), achieve simultaneously its channel and network 
capacities when there is at least one unstable mode of 
oscillation (time-scale) that is out of synchrony (see Sec. 

ixvm . 

While non-self-excitable channels permit the exchang- 
ing of a moderate amount of information in a reliable 
fashion, due to the low level of desynchronization in the 
channel, self-excitable channels permit the exchange of 
surprisingly large amounts of information, not necessary 
reliable, due to the higher level of desynchronization in 
the channel. 

In aiming at finding optimal network topologies, net- 
works that can not only transmit large amounts of in- 
formation but are also robust under alterations in the 
coupling strengths, we arrive at two relevant eigenvalues 
conditions which provide networks that satisfy all the 
optimal requirements. Either the network has elements 
that remain completely desynchronous for large varia- 
tions of the coupling strength, forming the self-excitable 
channels, or the network has elements almost completely 
synchronous, forming the non-self-excitable channels. In 
fact, the studied network, a network formed by electri- 
cally connected Hindmarsh-Rose neurons, can have si- 
multaneously self-excitable and non-self-excitable chan- 
nels. 

Self-excitable networks, namely those that have a ma- 
jority number of self-excitable channels, have the topol- 
ogy of a perturbed star, i.e., they are composed of a cen- 
tral neuron connected to most of the other outer neurons, 
and some outer neurons sparsely connected among them- 



selves. The networks that have non-self-excitable chan- 
nels have the topology of a perturbed fully connected 
network, i.e., a network whose elements are almost all- 
to-all connected. The self-excitable network has thus a 
topology which can be considered to be a model for mini- 
columnar structure of the mammalian neocortex [13| . 

In order to construct optimal networks, we have used 
two approaches. Firstly, (Sec. IVIIIAp . we use a 
Monte Carlo evolution technique [1J| to find the topol- 
ogy of the network, assuming equal bidirectional coupling 
strengths. This evolving technique simulates the rewiring 
of a neuron network that maximizes or minimizes some 
cost function, in this case a cost function which produces 
optimal networks to transmit information. In the sec- 
ond approach (Sec. IVIIIBp . we allow the elements to be 
connected with different coupling strengths. We then use 
the Spectral Theorem to calculate the coupling strengths 
of an all-to-all topology network. 

Finally, we discuss how to extend these results to net- 
works formed by elements that are non-chaotic (Sec. IIX[) . 
and to non-autonomous networks, that are being per- 
turbed by some time-dependent stimuli (Sees. IIXI and 

ig. 



IV. RESULTS 

V. UPPER BOUND FOR THE MUTUAL 

INFORMATION RATE (MIR) IN AN ACTIVE 

NETWORK 

In a recent publication [ll|, we have argued that the 
mutual information rate (MIR) between two elements in 
an active chaotic network, namely, the amount of infor- 
mation per unit time that can be realized in one element, 
fc, by measuring another element, /, regarded as Ic, is 
given by the sum of the conditional Lyapunov exponents 
associated with the synchronization manifold (regarded 
as A'') minus the positive conditional Lyapunov expo- 
nents associated with the transversal manifold (regarded 
as A^). So, /c = Ail -A-L. 

As shown in [l2|, if one has N=2 coupled chaotic sys- 
tems, which produce at most two positive Lyapunov ex- 
ponents Ai, A2 with Ai > A2, then A" = Ai and A-*^ — \2- 
Denote the trajectory of the element k in the network by 
Xk- For larger number of elements, N, the approaches 
proposed in [ll| remain valid whenever the coordinate 
transformation Xidii = Xk + xi (which defines the syn- 
chronization manifold) and Xki_L = Xk — xi (which de- 
fines the transversal manifold) successfully separates the 
two systems k and I from the whole network. Such a 
situation arises in networks of chaotic maps of the inter- 
val connected by a diffusively (also known as electrically 
or linear) all-to-all topology, where every element is con- 
nected to all the other elements. These approaches were 
also shown to be approximately valid for chaotic networks 
of oscillators connected by a diffusively all-to-all topol- 
ogy. The purpose of the present work is to extend these 



approaches and ideas to active networks with arbitrary 
topologies. 

Consider an active network formed by N equal ele- 
ments, Xi (i = 1,...,A'^), where every D-dimensional 
element has a different set of initial conditions, i.e., 
Xi 7^ X2 ^ . . . 7^ x^r . The network is described by 



F(x,)-a^a,,H(x,), 



(1) 



where Qij is the ij element of the coupling matrix. Since 
we choose ^ Gij = in order for a synchronization man- 
ifold to exist by the subspace 77 = Xi = X2 = X3 = . . . = 
xat, we can call this matrix the Laplacian matrix. The 
synchronous solution, 77, is described by 



77 = F{r,) 



(2) 



The way small perturbations propagate in the network 
|17| is described by the i {i — 1, . . . , N) variational equa- 
tions of Eqs. (HI), namely writing Xi = 77 -|- (5xi and ex- 
panding Eq. ([T]) in 5xi, 



S±, = [VF(x,) 



N 



(3) 



obtained by linearly expanding Eq. ^. The spectra of 
Lyapunov exponent is obtained from Eq. ([3]). 

Making x.^ — ^, which can be easily numerically done 
by setting the elements with equal initial conditions and 
taking H(xi) — Xi, Eq. ^ can be made block diagonal 
resulting in 



i^ = [VF(x,) - a7,]e» 



(4) 



where % are the eigenvalues (positive defined) of the 
Laplacian matrix ordered such that 'ji+i > 7^. Note that 
71=0. 

Notice that the network dynamics is described by Eq. 
dl]), which assumes that every element has different ini- 
tial conditions and therefore different trajectories (except 
when the elements are completely synchronized) . On the 
other hand, Eq. ^ that provides the conditional expo- 
nents considers that all the initial conditions are equal. 
The equations for ^1 describe the propagation of pertur- 
bations on the synchronization manifold ^, and the other 
equations describe propagation of perturbations on the 
manifolds transversal to the synchronization manifold. 
While Eq. ^ provides the set of Lyapunov exponents 
of an attractor, Eq. ^ provides the Lyapunov expo- 
nents of the synchronization manifold and its transversal 
directions. 

Notice also that when dealing with linear dynamics, 
the Lyapunov exponents [obtained from Eq. ([3])] are 
equal to the conditional exponents [obtained from Eq. 
(|3])] independently on the initial conditions. 

Then, the upper bound of the MIR that can be mea- 
sured from an element Xk by observing another element 



xi, i.e. the upper bound of the MIR in the communica- 
tion channel c'~^ is 



< lA^-A* 



(5) 



with i £ (2, . . . , N), and A* representing the sum of all 
the positive Lyapunov exponents of the equation for the 
mode £,i, in Eq. ([4]). So, A^ is the sum of the positive 
conditional exponents obtained from the separated varia- 
tional equations, using the smallest eigenvalue associated 
with the exponential divergence between nearby trajec- 
tories around ^, the synchronous state, and A* (z > 1) 
are the sum of the positive conditional exponents of one 
of the possible desynchronous oscillation modes. Each 
eigenvalue ji produces a set of conditional exponents A*„, 
with m = 1, . . . , _D. 

Each oscillatory mode ^^ represents a subnetwork 
within the whole network which possesses some oscilla- 
tory pattern. This oscillatory subnetwork can be used 
for communications purposes. Each mode represents a 
path along which information can be transmitted through 
the network. The oscillation mode associated with the 
synchronization manifold (^i) propagates some informa- 
tion signal everywhere within the network. The desyn- 
chronous modes limits the amount of information that 
one can measure from the signal propagated by the syn- 
chronous mode. Although Eq. ^ gives the upper bound 
for the amount of information between modes of oscil- 
lation, for some simple network geometries, as the ones 
studied here, we can relate the amount of information ex- 
changed between two vibrational modes to the amount 
of information between two elements of the network, and 
therefore, Eq. ([5]) can be used to calculate an upper 
bound for the MIR exchanged between pairs of elements 
in the network. For larger and complex networks, this 
association is non-trivial, and we rely on the reasonable 
argument that a pair of elements in an active network 
cannot transmit more information than some of the i — 1 
values of Pp^- 

The inequality in Eq. ([5]) can be interpreted in the fol- 
lowing way. The right hand side of Eq. ^ calculates 
the amount of information that one could transmit if 
the whole network were completely synchronous with the 
state ^, which is only true when complete synchroniza- 
tion takes place. Typically, we expect that the elements 
of the network will not be completely synchronous to ^. 
While the positive conditional exponents associated with 
the synchronization manifold provide an upper bound for 
the rate of information to be transmitted, the transversal 
conditional exponents provide a lower bound for the rate 
of erroneously information exchanged between nodes of 
the network. Thus, the amount of information provided 
by the right part of Eq. ([5]) overestimates the exact MIR 
which, due to desynchronization in the network, should 
be smaller than the calculated one. For more details on 
the derivations of Eq. ([5]), see Sec. IXVII 

Equation ^ allows one to calculate the MIR between 
oscillation modes of larger networks with arbitrary topol- 
ogy rescaling the MIR curve {Ip vs. a) obtained from two 



coupled elements. Denoting a*{N = 2) as the strength 
value for which the curve for A^ reaches a relevant value, 
say, its maximum value, then the coupling strength for 
which this same maximum is reached for A' in a network 
composed by N elements is given by 



'(N) 



2a* {N ^2) 



(6) 



where ^i{N) represents the ith largest eigenvalue of the 
A^-elements network. If the network has an all-to-all 
topology, thus, a* {N — 2) represents the strength value 
for which the curve of Ip reaches a relevant value, and 
a*{N) the strength value that this same value for Ip is 
reached. 

Notice that symmetries in the connecting network 
topology leads to the presence of degenerate eigenval- 
ues (=equal eigenvalues) in the Laplacian matrix, which 
means that there are less independent channels of com- 
munication along which information flows. Calling Q the 
number of degenerate eigenvalues of the Laplacian ma- 
trix, Eq. ^ will provide N — Q different values. 

As the coupling strength a is varied, the quantities 
that measure information change correspondingly. For 
practical reasons, it is important that we can link the 
way these quantities (see Sec. IXVp change with the 
way the different types of synchronization show up in 
the network (see Sec. IXVIII) . In short, there are three 
main types of synchronization observed in our examples: 
burst phase synchronization (BPS), when at least one 
pair of neurons are synchronous in the slow time-scale 
but desynchronous in the fast time-scale, phase synchro- 
nization (PS), when all pairs of neurons are phase syn- 
chronous, and complete synchronization (CS), when all 
pairs of neurons are completely synchronous. The cou- 
pling strength for which these synchronous phenomena 
appear are denoted by gbps, <^pSi and ucs (with no 
superscript index). 

Finally, there are a few more relevant coupling 
strengths, which characterize each communication chan- 
nel. First, cr^j„, for which the sum of the ith conditional 
exponents A' equals the value of A^. For a < crl^in^ the 
communication channel i (whose upper rate of informa- 
tion transmission depends on the two oscillation modes ^i 
and ^i) behaves in a self-excitable way, i.e., A^ < A*. For 
f > <^mim ^^ > A*. Secondly, (t** indicates the couphng 
strength at which PjT is maximal. Thirdly, a^g indi- 
cates the coupling strength for which the communication 



channel 



becomes "stable", i.e.. A* < 0. At ct = a" 



the self-excitable channel capacity of the channel c*^^ is 
reached and at ct = <^cs^ ^^^ non-self-excitable channel 
capacity is reached. Finally, ctc is the coupling for which 
the network capacity is reached, and then, when the KS- 
entropy of the network is maximal. For other quantities, 
see Sec. EVl 



VI. THE MIR IN NETWORKS OF COUPLED 
HINDMARSH-ROSE NEURONS 

We investigate how information is transmitted in self- 
excitable networks composed of N bidirectionally cou- 
pled Hindmarsh-Rose neurons |1|: 

ii = yi + ixl -xl - Zi + li + a'^ g^j (xj) 



Vt = 1 - 5a;j - yt 

Zi = ~rzi -f 4r(a;i -f 1.6) 



(7) 



The parameter r modulates the slow dynamics and is set 
equal to 0.005, such that each neuron is chaotic. The 
index i ^ j assumes values within the set [1 , . . . , iV] . 
Sk represents the subsystem formed by the variables 
ixk,yk, Zk) and Si represents the subsystem formed by 
the variables (xi,yi,zi), where fc=[l,...,A^ — 1] and 
l=[k + 1, . . . , A^]. The Laplacian matrix is symmetric, 
so Gji = Qij, and aQji is the strength of the electrical 
coupling between the neurons, and we take for /j the 
value li =3.25. 

In order to simulate the neuron network and to cal- 
culate the Lyapunov exponents through Eq. ([3]), we 
use the initial conditions x=-l.3078+ri, y— -7.3218+7], 
and 2=3.3530+77, where 77 is an uniform random number 
within [0,0.02]. To calculate the conditional Lyapunov 
exponents, we use the equal initial conditions, a;=-1.3078, 
2/=-7.3218, and 2=3.3530. 

All-to-all coupling: Here, we analyze the case where 
N neurons arc fully connected to every other neuron. The 
Laplacian matrix has A^ eigenvalues, 71 =0, and A^ — 1 de- 
generate ones "fi=N, i = 2, . . . , N. Every pair of neurons 
exchange an equal amount of MIR. Although, there are 
N X {N — l)/2 pairs of neurons, there is actually only 
one independent channel of communication, i.e., a per- 
turbation applied at some point of the network should be 
equally propagated to all other points in the network. In 
Fig. [TJA), we show the MIR, Ic, calculated using the ap- 
proaches in Refs. [ll|, [l3i ^Py calculated using the right 
hand-side of Eq. (O, and Is, calculated encoding the 
trajectory between pair of neurons (Sec. IXIIip . and the 
Kolmogorov-Sinai entropy, Hks, for a network composed 
by N=2 neurons. In (B), we show these same quantities 
for a network formed by N=A neurons. 

While for CT = and a > acs, we have that Ic ^ Ip = 
Is, for a = a^* (when the self-excitable channel capacity 
is reached) it is clear that Ip should be an upper bound 
for the MIR, since not only Ip > Ic but also Ip > Is- 
Notice the good agreement between Ic and Is, except 
for a = (T^,j„, when Is > Hks, which violates Eq. (fT3| . 

The star symbol indicates the value of the coupling, 
o'BPS (Sec. IXVII[) . for which burst phase synchronization 
(BPS) appears while the spikes are highly desynchronous. 
The appearance of BPS coincides with the moment where 
all the quantifiers for the MIR are large, and close to a 
coupling strength, ac, for which the network capacity is 
reached (when Hks is maximal). 



At this point, the network is sufficiently desynchronous 
to generate a large amount of entropy, which implies a 
large A*, for i > 2. This is an optimal configuration for 
the maximization of the MIR. There exists phase syn- 
chrony in the subspace of the slow time-scale z variables 
(which is responsible for the bursting-spiking behavior), 
but there is no synchrony in the (x, y) subspace. This 
supports the bindin g h ypothesis, a fundamental con- 
cept of neurobiology Q^ which sustains that neural net- 
works coding the same feature or object are functionally 
bounded. It also simultaneously supports the works of 
|15J . which show that desynchronization seems to play 
an important role in the perception of objects as well. 
Whenever A^ approaches zero, at cr = crcs, there is a 
drastic reduction in the value of Hks as well as Ip , since 
the network is in complete synchronization (CS), when 
all the variables of one neuron equals the variables of the 
other neurons. 

Therefore, for coupling strengths larger than the one 
indicated by the star symbol, and smaller than the one 
where CS takes place, there is still one time-scale, the 
fast time-scale, which is out of synchrony. 

For a > cr'^^^, the only independent communication 
channel is of the non-self-excitable type. That means 
A* < A^ (* > 2), and as the couphng strength increases, 
Hks decreases and Ip increases. 

Note that the curve for Ip shown in Fig. [U^B) can 
be obtained by rescaling the curve shown in Fig. HJA), 
applying Eq. ^. 

Nearest-neighbor coupling: Here, every neuron is 
connected to its nearest neighbors, with periodic bound- 
ary conditions, forming a closed ring. The eigenval- 
ues of the Laplacian matrix can be calculated from 

7fc=4sin(2l^^) , k G [l,...,N]. Notice that in this 
example, 7fc+i might be smaller than 7^ due to the de- 
generacies. We organize the eigenvalues in a crescent or- 
der. For our further examples, we consider N=A [in Fig. 
HA)] and 7V=6 [in Fig. ^B)]. For 7V=4, 71=0, 72.3=2, 
74=4, and for A^=6, 7i=0, 72,3=1, 74,5=3, 76=4. 

Networks with a nearest-neighbor coupling topology 
and an even number of elements possess a connecting ma- 
trix Q with N/2 — 1 degenerate eigenvalues, and therefore, 
N — N/2 + 1 distinct eigenvalues. There are only N — N/2 
different minimal path lengths connecting the elements 
of the network. The minimal path length quantifies the 
minimal distance between an element and another in the 
network by following a path formed by connected ele- 
ments. Note that Ip assumes only N — N/2 different 
values. It is reasonable to state that each different value 
corresponds to the exchange of information between ele- 
ments that have the same minimal path length. 

For a network with N—A [Fig. EI^A)], there are two pos- 
sible minimal path lengths, 1 and 2. Either the elements 
are 1 connection apart, or 2 connections apart. For such 
a network, it is reasonable to associate Ip = X^ — A^ with 
the MIR between two elements, Sk and Sk+2, that are 
2 connections apart, and Ip = |A^ ^ A^| to the MIR be- 
tween two elements, Sk and Sk+i, that are 1 connection 





FIG. 1; The quantities Ic (black circles), Ip (red squares), 
Is (green diamonds), and Hks (blue diamonds), for two (A) 
and four (B) coupled neurons, in an all-to-all topology. Notice 
that since there are only two different eigenvalues, there is 
only one channel of communication whose upper bound for 
the MIR is given by Ip = |A^ — A^|. Also, Is and Ic represent 
the mutual information exchanged between any two pairs of 
elements in the system. In (A), cr^* =0.092, cjbps — 0.2, 
o-Ln =0.42, aps = 0.47, and crcs=0.5. In (B), a^*=0.046, 
asps = 0.1, cr^„ =0.21, aps = 0.24, and acs=0.25. CS 
indicates the coupling interval a > acs for which there exists 
complete synchronization. 



apart. The more distant (closer) an element is from any 
other, the larger (smaller) the coupling strength for them 
to synchronize. In addition, a'^* > a'^* and cr^j„ > cr^i„. 
That means that the more distant elements are from each 
other the larger the coupling strength is, in order for 
these two elements to exchange a large rate of informa- 
tion, since a^* > cr'^*. In addition, since cr^j„ > cr^i„, 
the communication channel responsible for the exchange 
of information between closer elements (the channel c^) 
becomes non-self-excitable for a smaller value of the cou- 
pling strength than the strength necessary to turn the 
communication channel responsible for the exchange of 
information between distant elements (the channel c^) 
into a non-self-excitable channel. Since the level of desyn- 
chronization in a non-self-excitable channel is low, then, 
closer elements can exchange reliable information for 
smaller coupling strengths than the strength necessary 
for distant elements to exchange reliable information. 
Note that due to the 1 degenerated eigenvalue, Ip=Ip, 



a'^*, and a. 



= a3 . . 
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A similar analysis can be 



done for the network N=6, whose results are shown in 
Fig-HB). 

The KS entropy of the network, Hks, is also shown 



FIG. 2; The quantities Ip and Hks for nearest-neighbor 
networks with N=4: (A) and A^=6 (B). In (A), o-^*=0.09. 



=0.42, aL:„=0.21, aSs=0.25, aeps = 0.18, 



In (B), o-2*=o.l8, cr''*=0.061. 



^ mm ^.^ij, ^ rmn 

o"ps=0.462, and acs=0.5. 

aL„=0.84, <7Ln=0.27, a^s^O.SS, asps = 0.23, <7ps=0.78, 
and (Tcs = 1.0. The stars point to where BPS first appears. 
CS indicates the coupling interval a > acs for which there 
exists complete synchronization. 



in (B), o-i,„=0.84 and at 



mm 

easily derived from Eq 

^4 



in this figure. In (A), cr^j„=0.42 and (t^j„=0.21, and 



0.275, values that can be 
Note that the values of 
a = (J^j^n in (A) [and a — Omim ™ (B)] are close to the 
parameter for which BPS in the slow time-scale is first 
observed in these networks (indicated by the star symbol 
in Fig. HI), asps = 0.18 [in (A)] and obps ^ 0.23 [in 
(B)]. At a = a^,^ [in (A)] and a ^ a^,,, [in (B)], also 
the quantities Ip and Hks are large. 

Another important point to be emphasized in these 
networks is that Aa]^g^ = acs — <^mm > regarded as the 
non-self-excitable robustness parameter for the commu- 
nication channels c% with j=3 for the network with A'^=4 
[in (A)] and z = 5 for the network with A^=6 [in (B)] is 
large. This is a consequence of the fact that the normal- 
ized spectral distance (NED), (7^—72)/-/^ is also large, for 
either z = 4 [in (A)] or i = 6 [in (B)]. Having a large NED 
between the ith largest and the first largest eigenvalues 
results in a non-self-excitable channel, d~^ , robust under 
large alterations of the coupling strength. On the other 
hand, Acr^^ = cr^^jj, regarded as the self-excitable ro- 
bustness parameter for the communication channel c*~^, 
is large, for i = 2,3. This is a consequence of the fact that 
the normalized spectral distance (NED), (7JV — Ji)/N is 
large. Having a large NED between the largest and the 
ith largest eigenvalues results in a self-excitable chan- 
nel, c*~^, robust under large alterations of the coupling 
strength. 

Notice also that the maximal values of Ip for the all-to- 



all and nearest-neighbor networks topologies is the same 
(see Figs. [T]and[2]). This shows that the maximum of 
Ip does not depend on the number, N , of elements in 
the network. Not so in the case of the network capac- 
ity Cc^ which increases with N . Thus, pairs of elements 
can transmit information in a rather limited rate, but de- 
pending on the number of elements forming the network, 
a large number of channels can simultaneously transmit 
information. 

Star coupling: We consider iV=4. There is a cen- 
tral neuron, denoted by Si, bidirectionally connected to 
the other three {Sk,k = 2,3,4), but none of the others 
are connected among themselves. The eigenvalues of the 
Laplacian matrix are 71=0,72,3=1,74 = N. So, not only 
the NED between jn and jn-i is large but also between 



7Ar and Jn-2, and therefore, Aa^^ and Aa^^ are 
large. This provides a network whose channels c^ and c^ 
have a large MIR for a large coupling strength alteration. 
Note that if 7Ar_i is far away from jn that implies that 
7Ar_2 is also far away from 7jv. Thus, a reasonable spec- 
tral distance between jn-i and jn is a "biological re- 
quirement" for the proper function of the network, since 
even for larger coupling strengths there will be at least 
one oscillation mode which is desynchronous, a configura- 
tion that enables perturbation (meaning external stimuli) 
to be propagated within the network [l6|. 
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FIG. 3: MIR between the central neuron and an outer one 
(black circles), Ip, (resp. /s(l, k), in green line), and between 
two outer ones (red squares), Ip, (resp. Is{k, I), in blue line). 
Blue diamonds represents the KS-entropy. Other quantities 
are a"* = 0.181, a^* = 0.044, at.in = 0.84, o-^n = 0.22, 
0-^3=0-27, o-i3Ps=0.265, crps=0.92, and crcs=l-0. The star 
indicates the parameter for which BPS first appears. 

The largest eigenvalue is related to an oscillation mode 
where all the outer neurons are in synchrony with each 
other but desynchronous with the central neuron. So, 



here it is clear the association between jA"'^ — A^| and the 
MIR between the central neuron with an outer neuron, 
since A^ represents the amount of information of the syn- 
chronous trajectories among all the neurons, while A'* is 
the amount of information of the desynchronous trajec- 
tories between the central neuron and any outer neuron. 
The other eigenvalues (72,73) represent directions trans- 
verse to the synchronization manifold in which the outer 
neurons become desynchronous with the central neuron 
in waves wrapping commensurately around the central 
neuron [l3|. Thus, A^ and A'^ are related to the error 
in the transmission between two outer neurons, k and I, 
with k,l ^ 1. 

Note that the MIR between Si and an outer neuron 
(upper bound represented by /|, = |A^ — A"*] and Is rep- 
resented by Is{l,k), in Fig. [3]) is larger (smaller) than 
the MIR between two outer neurons (upper bound rep- 
resented by /p = |A^ — A^l and Is is represented by 
Is{k,l), in Fig. [3]), for small coupling (for when the 
channel C^ is self-excitable, and a > cr^jn). Similar to 
the nearest-neighbor networks, the self-excitable and the 
non-self-excitable channel capacities of the channel as- 
sociated with the transmission of information between 
closer elements (the channel c^) are achieved for a smaller 
value of the coupling strength than the one necessary to 
make the channels associated with the transmission of 
information between more distant elements (the channel 
c^) to achieve its two channel capacities. That property 
permits this network, for a = '^mim ^'^ transmits simul- 
taneously reliable information using the channel c^ and 
with a higher rate using the channel c^ . 

Notice, in Fig. [31 that a'^* = af^^„ ~ ctbps — o^c- So, 
when the channel capacity of the channel c^ is reached, 
also Hxs of the network is maximal, and the network 
operates with its capacity. 

Another point that we want to emphasize in this net- 
work is that while a large NED between 7jv and jn-i 
provides a network whose channel c^ is self-excitable and 
can transmit information at a large rate for a large cou- 
pling strength interval, a large NED between 73 and 72 
leads to a non-self-excitable channel c"^ even for small val- 
ues of the coupling amplitudes, and it remains non-self- 
excitable for a large variation of the coupling strength. 
Thus, while a large NED between the second and the first 
largest eigenvalues leads to a network whose channels are 
predominantly of the self-excitable types, a large NED 
between the second largest and the third largest eigen- 
values provide a network whose communication channels 
are predominantly of the non-self-excitable types. 



VII. EIGENVALUES CONDITIONS FOR 
OPTIMAL NETWORK TOPOLOGIES 

Finding network topologies and coupling strengths in 
order to have a network that operates in a desired fashion 
is not a trivial task (see Sec. IXVIIII and IXIXp . An ideal 
way to proceed would be to evolve the network topology 



in order to achieve some desired behavior. In this paper, 
we are interested in maximizing simultaneously /p, the 
KS-entropy, and the average (Ip), for a large range of 
the coupling strength, characteristics of an optimal net- 
work. However, evolving a network in order to find an 
optimal one would require the calculation of the MIR in 
every communication channel and Hks for every evo- 
lution step. For a typical evolution, which requires 10^ 
evolution steps, such an approach is impractical. 

Based on our previous discussions, however, an optimal 
network topology can be realized by only selecting an 
appropriate set of eigenvalues which have some specific 
NED. Evolving a network by the methods of Sees. IXVIIII 



SE 



and 1X1X1 using a cost function which is a function of only 
the eigenvalues of the Laplacian matrix is a practical and 
physible task. 

The present section is dedicated to describe the deriva- 
tion of this cost function. 

We can think of two most relevant sets of eigenvalues 
which create optimal networks, and they are represented 
in Fig. Ul Either it is desired eigenvalues that produce 
a network predominantly self-excitable [SE, in Fig. 2] or 
predominantly non-self-excitable [NSE, in Fig. 0]. 

In a network whose communication channels are pre- 
dominantly self-excitable, it is required that the NED 
ilN — ^N-i)/N is maximal and {■yN-i)/N minimal. 
Therefore, we want a network for which the cost func- 
tion 



Bi 



In — lN~i 

IN-I 



(8) 



is maximal. 

A network whose eigenvalues maximize Bi has self- 
excitable channels for a large variation of the coupling 
strength. As a consequence, {Ip) as well as H^s is large 
for ere [cr^„,cr^ J. 

In a network whose communication channels are pre- 
dominantly non-self-excitable, it is required that the 
NED (73 — 72)/-/V is maximal and {"f2)/N minimal. 
Therefore, we want a network for which the cost func- 
tion 



B2 



73 -72 
72 



(9) 



is maximal. 

A network whose eigenvalues maximize the condition 
in Eq. ([9]) have non-self-excitable channels for a large 
variation of the coupling strength. As a consequence, 
(Ip) is large for a € [o'minj'^minli which is a small cou- 
pling range, but since there is still one oscillation mode 
that is unstable (the mode ^2), Hks is still large for a 
large range of the couphng strength (cr < a'^in)- Most 
of the channels will transmit information in a reliable 
way, since the error in the transmission, provided by A* 
{i > 2), of most of the channels will be zero, once A* < 0. 

Since degenerate eigenvalues produce networks with 
less vibrational modes, we assume in the following the 
absence of such degenerate eigenvalues. In addition, we 
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FIG. 4: Representation of the eigenvalues sets that produce 
optimal self-excitable (SE) and non-self-excitable active net- 
works (NSE). 



assume that there is a finite distance between eigenvalues 
so that the network becomes robust under rewiring, and 
therefore, perturbing Qij will not easily create degenerate 
eigenvalues. 

A network that is completely synchronous and has no 
unstable modes does not provide an appropriate envi- 
ronment for the transmission of information about an 
external stimulus, because they prevent the propagation 
of perturbations. Networks that can be easily completely 
synchronized (for small coupling strengths) requires the 
minimization of 77V— 72, or in terms of the eigenratio, the 
minimization of "fN/"f2- We are not interested in such a 
case. To construct network topologies that are good for 
complete synchronization, see Refs. jlTI. llSl. [l9j . 



VIII. OPTIMAL TOPOLOGIES FOR 
INFORMATION TRANSMISSION 

Before explaining how we obtain optimal network 
topologies for information transmission, it is important 
to discuss the type of topology expected to be found by 
maximizing either Bi, in Eq. ([S]) or S2, in Eq. Q. No- 
tice that Laplacians whose eigenvalues maximize Bi are 
a perturbed version of the star topology, and the ones 
that maximize B2 are a perturbed version of the all-to- 
all topology. In addition, in order to have a network that 
presents many independent modes of oscillations it is re- 
quired that the Laplacian matrix presents as much as 
possible, a large number of non-degenerate eigenvalues. 
That can be arranged by rewiring (perturbing) networks 
possessing either the star or the nearest-neighbor topol- 
ogy, breaking the symmetry. 



In order to calculate an optimal Laplacian, we propose 
two approaches. 

One approach, described in Sec. IXVIIIl is based on 
the reconstruction of the network by evolving techniques, 
simulating the process responsible for the growing or 
rewiring of real biological networks, a process which tries 
to maximize or minimize some cost function. The results 
are discussed in Sec. IVIII Al 

A second approach, described in Sec. IXIXl is based 
on the Spectral Theorem, and produces a network in or- 
der for its Laplacian to have a previously chosen set of 
eigenvalues. These eigenvalues can be chosen in order to 
maximize the cost function. The results are discussed in 
Sec. IVIII Bl 



A. Evolving networks 

In order to better understand how a network evolves 
(grows) in accordance with the maximization of the cost 
functions in Eqs. ^ and ^, we first find the network 
configurations with a small number of elements. To be 
specific, we choose iV=8 elements. To show that in- 
deed the calculated network topologies produce active 
networks that operate as desired, we calculate the aver- 
age upper bound value of the MIR [Eq. (|12p] for neu- 
ral networks described by Eqs. ([7]) with the topology 
obtained by the evolution technique, and compare with 
other network topologies. Figure [S] shows (/p), the aver- 
age channel capacity, calculated for networks composed 
of 8 elements, using one of the many topologies obtained 
by evolving the network maximizing Bi (circles, denoted 
in Fig. by "evolving 1"), all-to-all topology (squares), 
star topology (diamonds), nearest-neighbor (upper trian- 
gle), and maximizing B2 (down triangle, denoted in Fig. 
by " evolving 2" ) . The star points to the value of (T^j„ , 
when c^, the most unstable communication channel (a 
self-excitable channel), becomes non-self-excitable. 

As desired the evolving network 1 has a large upper 
bound for the MIR (as measured by (Ip)) for a large 
range of the coupling strength, since the network has 
predominantly self-excitable channels. The channel c^ 
has a large robustness parameter Acr|g, i.e., it is a self- 
excitable channel for a < cr,^j„, where cr^i„=2.0. In 
contrast to the other topologies, in the star, nearest- 
neighbor, and all-to-all topologies, A(t|^ is smaller and 
Aaffg^ is larger. Even though most of the channels in 
the evolving 2 topology are of the non-self-excitable type, 
(Ip) remains large even for higher values of the coupling 
strength. That is due to the channel c^ which turns into 
a self-excitable channel only for a > 2. 

The KS-entropies of the 5 active networks whose (Ip) 
are shown in Fig. [5] are shown in Fig. [51 Typically, the 
network capacities are reached for roughly the same cou- 
pling strength for which the maximum of (Ip), is reached. 
In between the coupling strength for which the network 
capacities and the maximal of {Ip) are reached, A'^ be- 
comes negative. At this point, also BPS appears in the 
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FIG. 5: The average value of the upper bound MIR, (Ip) 
[as defined in Eq. (I12|l ] for active networks composed of 8 
elements using one of the many topologies obtained by evolv- 
ing the network maximizing Si (circles), all-to-all topology 
(squares), star topology (diamonds), nearest- neighbor (up- 
per triangle), and maximizing B2 (down triangle). The val- 
ues of a^in indicated by the starts are (t^j„=0.169 (evolv- 



ing 1), 



=0.05 (all-to-aU), a4i„=0.037 (star), a^i„=0.037 



(nearest-neighbor), and (T^i„=0.6 (evolving 2). The evolv- 
ing 1 network has a Laplacian with relevant eigenvalues 
77=3.0000, 78=6.1004, which produces a cost function equal 
to Bl =1.033. The evolving 2 network has a Laplacian with rel- 
evant eigenvalues 72=0.2243 and 73 = 1.4107, which produces 
a cost function equal to 02=5.2893. 



slow time-scale, suggesting that this phenomena is the 
behavioral signature of a network that is able to trans- 
mit not only large amounts of information between pairs 
of elements (high MIR) but also overall within the net- 
work (high Hks)- 

Note however, that since the evolving networks have 
a small number of elements, the cost function cannot 
reach higher values and therefore, the networks are not 
as optimal as they can be. For that reason, we proceed 
now to evolve larger networks, with 7V=32. 

Maximization of the cost function Bi leads to the net- 
work connectivity shown in Fig. [7KA) and maximization 
of the cost function B2 leads to the network connectivity 
shown in Fig. [TJB). In (A), the network has the topol- 
ogy of a perturbed star, a neuron connected to all the 
other outer neurons, thus a hub, and each outer neuron 
is sparsely connected to other outer neurons. The arrow 
points to the hub. In (B),the network has the topology of 
a perturbed all-to-all network, where elements are almost 
all-to-all connected. Note that there is one element, the 
neuron S'32, which is only connected to one neuron, the 
Si. This isolated neuron is responsible to produce the 
large spectral gap between the eigenvalues 73 and 72. 
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FIG. 6; KS-entropy for the same active networks of Fig. [5] 
composed of 8 elements. 



(/p) for the network topology represented in Fig. [Tl^A) 
is shovifn in Fig. [8] as circles, and (/p) for the net- 
work topology represented in Fig. EI^B) is shown in 
Fig. [5] as squares. We see that the star topology, 
whose connectivity is represented in [ZIA), has larger 
(/p) for a larger coupling strength than the topol- 
ogy whose connectivity is represented in [7]JB). Other 
relevant parameters of the network whose topology is 
represented in [TjA) are cr^j„=0.8468, o-^i„=0.8249, 
, cri!,„ =0.0278, crcs=0.9762 and for the topology 

=0.8512, aL„=0.042, 



2 

rain 



represented in [7{^B) are a: 
cr^^^=0.031, and crcs=0.9761. 

It is worth to comment that the neocortex is being 
simulated in the Blue Brain project, by roughly creat- 
ing a large network composed of many small networks 
possessing the star topology. By doing that, one tries to 
recreate the way minicolumnar structures |13j are con- 
nected to minicolumnar structures of the neocortex [20[ . 
Each minicolumn can be idealized as formed by a pyra- 
midal neuron (the hub) connected to its interneurons, 
the outer neurons in the star topology, which are respon- 
sible for the connections among this minicolumn (small 
network) to others minicolumn. So, the used topology 
to simulate minicolumns is an optimal topology in what 
concerns the transmission of information. 



B. Constructing a network by a given set of 
eigenvalues 

It is of general interest to assess if the eigenvalues 
obtained from the method in Sec. IXVIIII (in order to 
have a network Laplacian whose eigenvalues maximize 
the cost function B) can be used to construct other net- 



FIG. 7: A point in this figure in the coordinate k X I means that 
the elements 5fe and Si are connected with equal couplings in a 
bidirectional fashion. In (A), a 32 elements network, constructed 
by maximizing the cost function S\ in Eq. JSjl and in (B), 32 el- 
ements network, constructed by maximizing the cost function B2 
in Eq. l(9]l. In (A), the network has the topology of a perturbed 
star, a hub of neurons connected to all the other neurons, where 
each outer neuron is sparsely connected to other neurons. The ar- 
row points to the hub. In (B),the network has the topology of a 
perturbed all-to-all network, where elements are almost all-to-all 
connected. Note that there is one element, the neuron ^32, which 
is only connected to one neuron, the Si. This isolated neuron is 
responsible to produce the large spectral gap between the eigenval- 
ues 73 and 72. In (A), the relevant eigenvalues are 731=4.97272, 
732=32, which produce a cost function equal to 6=5.43478. In 
(B), the relevant eigenvalues are 72=0.99761, 73=27.09788, which 
produce a cost function equal to 02=26.1628. 



works (whose Laplacian preserve the eigenvalues) main- 
taining still the properties here considered to be vital for 
information transmission. 

By a given set of eigenvalues, one can create a Lapla- 
cian matrix, Q' ^ with non-zero real entries, using the 
method described in Sec. IXIXI The resulting network 
will preserve the eigenvalues and the synchronous solu- 
tion in Eq. ([2]) , which means that the values of /p of the 
topology created by the method in Sec. IXIXI are equal 
to the values of the network topologies that provide the 
set of eigenvalues, in the following example, the network 
connectivities represented in Fig. [ZKA-B). 

In Fig. IHA-B), circles represent the values of Hks for 
the network whose connectivities are represented in Figs. 
[71[A-B), and the squares represent this same quantity for 
a network whose Laplacian is calculated by the method 
in Sec. IXIXI in order to preserve the same eigenvalues of 
the network topologies represented in Fig. EJA-B). 

The main difference between these two networks is that 
for the one constructed by the method in Sec. IXIXI when 
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FIG. 8: {Ip) for the networks shown in Fig. EJA-B) by circles 
and squares, respectively. 




FIG. 9: Hks for a network evolved by the method in Sec. 
IXVIIII in circles, and in squares, for a network whose Lapla- 
cian is calculated by the method in Sec. IXIXl in order for the 
Laplacian to generate the same eigenvalues as the ones gen- 
erated by the network Laplacian calculated by the evolution 
technique. In (A), we consider the same network topology 
whose connectivity is represented in Fig. [ZlA), and in (B), 
we consider the same network topology whose connectivity is 
represented in Fig. E^B). Note that in general, a Laplacian of 
an active network whose elements are connected with different 
coupling strengths, possess a smaller value of Hks- 



X^ becomes zero, simultaneously A^ becomes also zero, a 



consequence of the fact that all the neurons enter in a 
non-trivial but periodic oscillation. In general, however, 
both networks preserve the characteristics needed for op- 
timal information transmission: large amounts of MIR 
and HkSi however, the ones constructed by the evolution 
technique have larger HkSi ^nd possess a larger MIR for 
larger ranges of the coupling strength. The network ob- 
tained by the method in Sec. IXIXl is more synchronizable, 
a consequence of the fact that the coupling strengths are 
non-equal [ll[li]. 



IX. ACTIVE NETWORKS FORMED BY 

NON-CHAOTIC ELEMENTS 

The purpose of the present work is to describe how in- 
formation is transmitted via an active media, a network 
formed by dynamical systems. There are three possi- 
ble asymptotic stable behaviors for an autonomous dy- 
namical system: chaotic, periodic, or quasi-periodic. A 
quasi-periodic behavior can be usually replaced by either 
a chaotic or a periodic one, by an arbitrary perturba- 
tion. For that reason, we neglect such a state and focus 
the attention on active channels that are either chaotic 
or periodic. 

The purpose of the present section is dedicated to an- 
alyze how a source of information can be transmitted 
through active channels that are non-chaotic, that is pe- 
riodic, and that possess negative Lyapunov exponents. 

Equation (O is defined for positive exponents. How- 
ever, such an equation can also be used to calculate 
an upper bound for the rate of mutual information in 
systems that also possess negative Lyapunov exponents. 
Consider first a one-dimensional contracting system be- 
ing perturbed by a random stimulus. Further consider 
that the stimulus changes the intrinsic dynamics of this 
system. This mimics the process under which an active 
element adapts to the presence of a stimulus. 

Suppose the stimulus, 0„, can be described by a dis- 
crete binary random source with equal probabilities of 
generating '0' or '1'. Whenever 6'„ = 0, the system 
presents the dynamics a;„_|_i = a;„/2, otherwise a;„-)-i = 
(1 + Xn)/2. It is easy to see that the only Lyapunov 
exponent of this mapping, Ai , which is equal to the con- 
ditional exponent, A^, is negative. Negative exponents do 
not contribute to the production of information. From 
Eq. ([5]) one would arrive at Ip—Q. However, all the 
information about the stimulus is contained in the tra- 
jectory. If one measures the trajectory a;„, one knows 
exactly what the stimulus was, either a '0' or a T'. The 
amount of information contained in the stimulus is log (2) 
per iteration which equals the absolute value of the Lya- 
punov exponent, |Ai|. In fact, it is easy to show that 
Ic = Ip — lA"*"! = |Ai| = log (2), or if we use the inter- 
pretation of [21|, Ic — Ip = A, where A = |Ai| is the 
positive Lyapunov exponent of the time-inverse chaotic 
trajectory, Xn+m,Xn+m-i, ■ ■ ■ ,Xo, which equals the rate 
of information production of the random source. So, in 
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this type of active communication channel, one would 
consider in Eq. ([5]) the positive Lyapunov exponents of 
the time-inverse trajectory, or the absolute value for the 
negative Lyapunov exponent. 

Another example was given in |12l |. In this reference 
we have shown that a chaotic stimulus perturbing an ac- 
tive system with a space contracting dynamics (a neg- 
ative Lyapunov exponent) might produce a fractal set. 
We assume that one wants to obtain information about 
the stimulus by observing the fractal set. The rate of 
information retrieved about the stimulus on this fractal 
set equals the rate of information produced by the fractal 
set. This amount is given by -Di | A| , where Di is the infor- 
mation dimension of the fractal set and |A| the absolute 
value of the negative Lyapunov exponent. In fact, -Di|A| 
is also the rate of information produced by the stimulus. 
So, if an active system has a space contracting dynam- 
ics, the channel capacity equals the rate of information 
produced by the stimulus. In other words, the amount 
of information that the system allows to be transmit- 
ted equals the amount of information produced by the 
chaotic stimulus. 



X. THE ROLE OF A TIME-DEPENDENT 
STIMULUS IN AN ACTIVE NETWORK 

The most general way of modeling the action of an ar- 
bitrary stimulus perturbing an active network is by stim- 
ulating it using uncorrelated white noise. Let us assume 
that we have a large network with all the channels oper- 
ating in non-self-excitable fashion. We also assume that 
all the transversal eigenmodes of oscillations except one 
are stable, and therefore do not suffer the influence of 
the noise. Let us also assume that the noise is acting 
only on one structurally stable {— far from bifurcation 
points) element, Sk- To calculate the upper bound of the 
MIR between the element Sk and another element Si in 
the network, we assume that the action of the noise does 
not alter the value of A^. Then, the noise on the element 
Sk is propagated along the vibrational mode associated 
with the one unstable transversal direction, whose con- 
ditional exponent is A^. As a consequence, the action 
of the noise might only increase A^, while not affecting 
the negativeness of all the other exponents (A™, to > 2), 
associated with stable transversal modes of oscillation. 
That means that the channels responsible for transmit- 
ing large amounts of information (associated with A™, 
with m large) will not be affected. So, for such types of 
noises, Eq. ^ of the autonomous network is an upper 
bound for the non- autonomous network. 

Consider now a situation where the noise acts equally 
on all the elements of an active network. The mapping 
of Eq. pij) was proposed as a way to understand such 
a case. Consider the non-self-excitable map for s=-l. 
Note that the term p{x^ + y"^) that enters equally in all 
the maps has statistical properties of an uniformly dis- 
tributed random noise. Calculating Ip for p = (the 



noise-free map) we arrive at Ip ~ 2(7, for small cr, while 
the true MIR Ic = 2{a — p). These results are con- 
firmed by exact numerical calculation of the Lyapunov 
exponents of Eq. p^ as well as the calculation of the 
conditional exponents of the variational equations. So, 
this example suggests that Eq. ([5]) calculated for an au- 
tonomous non-perturbed network gives the upper bound 
for the mutual information rate in a non-autonomous net- 
work. 



XI. DISCUSSIONS 

We have shown how to relate in an active network the 
rate of information that can be transmitted from one 
point to another, regarded as mutual information rate 
(MIR), the synchronization level among elements, and 
the connecting topology of the network. By active net- 
work, we mean a network formed by elements that have 
some intrinsic dynamics and can be described by classical 
dynamical systems, such as chaotic oscillators, neurons, 
phase oscillators, and so on. 

Our main concern is to suggest how to construct an op- 
timal network. A network that simultaneously transmits 
information at a large rate, is robust under couplings 
alterations, and further, it possesses a large number of 
independent channels of communication, pathways along 
which information travels. 

We have studied two relevant conditions that the eigen- 
values of the Laplacian matrix have to satisfy in order for 
one to have an optimal network. The Laplacian matrix 
describes the coupling strengths among each element in 
the active network. 

The two eigenvalues conditions are designed in order 
to produce networks that are either self-excitable [max- 
imizing Eq. ([5])] or non-self-excitable [maximizing Eq. 
di])] (see definition of self-excitability in Sec. IXIVp . Self- 
excitable networks have communication channels that 
transmit information in a higher rate for a large range 
of the coupling strength. Most of the oscillation modes 
in these networks are unstable, and therefore, informa- 
tion is mainly propagated in a desynchronous environ- 
ment. Non-self-excitable networks have communication 
channels that transmit information in a higher rate for a 
small range of the coupling strength, however, they have 
channels that transmit reliable information in a moder- 
ate rate for large range of coupling strengths. Most of the 
oscillation modes in these networks are stable, and there- 
fore, information is mainly propagated in a synchronous 
environment, a highly reliable environment for informa- 
tion transmission. 

Therefore, to determine the topology of an optimal 
network one does not need to know information about the 
intrinsic dynamics of the elements forming the network. 

Once the network topology is obtained such that the 
eigenvalues of the Laplacian matrix maximizes either the 
cost function in Eq. ^ or the one in Eq. ^ , the actual 
amount of information that can be transmitted using the 
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obtained topology will depend on the intrinsic dynamics 
of the elements forming the network [F in Eq. ([1])] and 
also on the type of coupling [H in Eq. ([T])], of only two 
coupled elements [see Eq. ^]. 

In the examples studied here, phase synchronization 
(PS) in the subspace {x,y) results in a great decrease 
of the KS-entropy (See Figs. [1] and [2]) as well as of the 
MIR and Ip. However, a special type of partial phase 
synchronization, the BPS, appears simultaneously when 
some communication channel achieves its capacity. So, 
BPS ;12] can provide an ideal environment for informa- 
tion transmission, probably a necessary requirement in 
the brain [23, [2^ . Similarly, in networks of Rossler os- 
cillators, a type of non-self-excitable network, PS is the 
phenomenon! responsible to identify when the network is 
operating in a regime of high MIR [Il|, [2J] . 

In order to construct an optimal network, we have used 
two approaches. One based on a Monte Carlo evolving 
technique, which randomly mutates the network topolo- 
gies in order to maximize the cost functions in Eqs. ([8]) 
and ([9]) (see Sec. IXVIII|) . We do not permit the existence 
of degenerate eigenvalues. As a consequence 7jv — 7Ar-i 
as well as 73—72 is never zero. The mutation is performed 
in order to maximize the cost function, but we only con- 
sider network topologies for which the value of the cost 
functions Bi and B2 remain constant for about 10,000 
iterations of the evolving technique, within 1,000,000 it- 
erations. Even though more mutations could lead to 
networks that have larger values of the cost function, 
we consider that a reasonably low number of mutations 
would recreate what usually happens in real networks. 
The other approach creates an arbitrary Laplacian which 
reproduces a desired set of eigenvalues. 

Although both topologies provide larger amounts of 
MIR and Has, meaning large network and channel ca- 
pacities, the topology provided by the evolution tech- 
nique, which consider coupling strengths with equal 
strengths, is superior in what concerns information trans- 
mission. That agrees with the results of Ref. jl2] which 
say that networks composed by elements with non-equal 
control parameters can transmit less information than 
networks formed by equal elements, since networks whose 
coupling strengths are non-equal can be considered to be 
a model for networks with non-equal control parameters. 

So, if brain-networks somehow grow in order to maxi- 
mize the amount of information transmission, simultane- 
ously remaining very robust under coupling alterations, 
the minimal topology that small neural networks must 
have should be similar to the one in Fig. [^A), i.e., a net- 
work with a star topology, presenting a central element, a 
hub, very well connected to other outer elements, which 
are sparsely connected. 

Even though most of the examples worked out here 
concern simulations performed in a neural network of 
electrically coupled Hindmarsh-Rose neurons, our theo- 
retical approaches to find optimal topologies can be used 
to a large class of dynamical systems, in particular also 
to networks of synaptically (chemically) connected neu- 



rons. A neural network with neurons connected chem- 
ically would also be optimal if one connect neurons by 
maximizing either Eq. ([5]) or Eq. ^. The novelty 
introduced by the chemical synapses is that it can en- 
hance (as compared with the electrical synapses) both 
the self-excitable (using excitable synapses) or the non- 
self-excitable (using inibitory synapses) characteristic of 
the communication channels as well as it can enhance 
{Ip) [23 • From the biological point-of-view, of course, 
the chemical synapses provide the long-range coupling 
between the neurons. So, the simulations performed here 
for the larger HR networks should be interpreted as to 
simulations of a general active network, since neurons 
connected electrically can only make nearest-neighbor 
connections. 



XII. METHODS 



XIII. CALCULATING THE MIR BY SYMBOLIC 
ENCODING THE TRAJECTORY 



The MIR between two neurons can be roughly esti- 
mated by symbolizing the neurons trajectory and then 
measuring the mutual information from the Shannon en- 
tropy [23] of the symbolic sequences. From [23|, the mu- 
tual information between two signals Sk and Si is given 
by 



I's 



HiSk)-HiSi\Sk) 



(10) 



H{Sk) is the uncertainty about what Sk has sent (en- 
tropy of the message), and H{Si\Sk) is the uncertainty 
of what was sent, after observing Si. In order to esti- 
mate the mutual information between two chaotic neu- 
rons by the symbolic ways, we have to proceed with a 
non-trivial technique to encode the trajectory, which con- 
stitutes a disadvantage of such technique to chaotic sys- 
tems. We represent the time at which the n-th spike 
happens in Sk by T^, and the time interval between the 
n-th and the (n-l-l)-th spikes, by ST^. A spike happens 
when Xk becomes positive and we consider about 20000 
spikes. We encode the spiking events using the follow- 
ing rule. The i-th symbol of the encoding is a "1" if 
a spike is found in the time interval [iA, (i + 1)A[, and 
"0" otherwise. We choose A e [min ((JT^'), max((5T^")] 
in order to maximize Ig. Each neuron produces a sym- 
bolic sequence that is split into small non-overlapping 
sequences of length L—12. The Shannon entropy of the 
encoding symbolic sequence (in units of bits), is esti- 
mated by max 77 1 = -X^m ^™ ^0S2 -fm where Pm is the 
probability of finding one of the 2^ possible symbolic se- 
quences of length L. The term II{Si\Sk) is calculated by 
H{Si\Sk)=-HiSi)+H{Sk; Si), with HiSk-, Si) represent- 
ing the Joint Entropy between both symbolic sequences 
for Sk and Si. 

Finally, the MIR (in units of bits/unit time), Is, is 
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calculated from 



Is^ 



I's 



Ax L 



(11) 



The calculation of the Is by means of Eq. (fTTj) should 
be expected to underestimate the real value for the MIR. 
Since the HR neurons have two time-scales, a large se- 
quence of sequential zeros in the encoding symbolic se- 
quence should be expected to be found between two 
bursts of spikes (large 5TJ} values) , which lead to a reduc- 
tion in the value of H{Sk) followed by an increase in the 
value of H{Si\Sk), since there will be a large sequence of 
zeros happening simultaneously in the encoding sequence 
for the interspike times of Sk and Si . 



XIV. SELF-EXCITABILITY 



the maximal value of the average value 



N 



1 



(^-)-Ew^\>^'-^ 



1=2 



fl2) 



The network capacity of a network composed of N el- 
ements, Cn{N), is defined to be the maximum value of 
the Kolmogorov-Sinai (KS) entropy, Hks, of the net- 
work. For chaotic networks, the KS-entropy, as shown 
by Pesin [28|, is the sum of all the positive Lyapunov 
exponents. Notice that if / denotes the MIR then 



I <H 



KS 



(13) 



As shown in Ref. [I4I and from the many examples 
treated here, Cn{N) oc A'^, and so, the network capacity 
grows linearly with the number of elements in an active 
network. 



In Ref. [12] self-excitability was defined in the follow- 
ing way. An active network formed by N elements, is 
said to be self-excitable if Hks{N, o") > Hks{N, ct = 0), 
which means that the KS-entropy of the network in- 
creases as the coupling strength is increased. Thus, for 
non self-excitable systems, an increase in the coupling 
strength among the elements forming the network leads 
to a decrease in the KS-entropy of the network. 

Here, we adopt also a more flexible definition, in terms 
of the properties of each communication channel. We 
define that a communication channel c* behaves in a self- 
excitable fashion if A* > A^. It behaves in a non-self- 
excitable fashion if A* < A^. 



XVI. UNDERSTANDING EQ. (HJ 

Let us study Eq. ([5]) using an analytical example. For 
an introduction to the quantities shown here see Sec. IXVI 
Consider the following two coupled maps: 



Xn+i = 2xn - pa;„ -I- 2sa{y„ 
yn+i = 2y„ - pyl + 2scr(x„ 



Vn), 



(14) 



with p > 0, s = ±1, and Xn,yn G [0,1], which can be 
accomplished by applying the mod{\) operation. 



XV. MUTUAL INFORMATION RATE (MIR), 

CHANNEL CAPACITY, AND NETWORK 

CAPACITY 

In this work, the rate with which information is ex- 
changed between two elements of the network is calcu- 
lated by different ways. Using the approaches of Refs. 
[111 . [12 ] , we can have an estimate of the real value of the 
MIR, and we refer to this estimate as Ic- Whenever we 
use Eq. ([5]) to calculate the upper bound for the MIR, 
we will refer to it as Ip. Finally, whenever we calculate 
the MIR through the symbolic encoding of the trajectory 
as described in Sec. IXIIIl we refer to it as Is- 

We define the channel capacity of a communication 
channel formed by two oscillation modes depending on 
whether the channel behaves in a self-excitable fashion or 
not. So, for the studied network, every communication 
channel possess two channel capacities, the self-excitable 
capacity and the non-self-excitable one. A channel d 
operates with its self-excitable capacity when Ip is max- 
imal, what happens at the parameter (t(*+^)*. It operates 
with its non-self-excitable capacity when A*+^ = 0. 

We also define the channel capacity in an average sense. 
In that case, the averaged channel capacity is given by 



A. Positiveness of the MIR in Eq. (fTJ]) 

Here, we assume that p=Q. This map produces two 
Lyapunov exponents Ai=log(2) and A2=log (2 — 4so-). 
Since this map is linear, the conditional exponents are 
equal to the Lyapunov exponents. 

Using the same ideas of Ref. [I4] , actually an interpre- 
tation of the way Shannon '27] defines mutual informa- 
tion, the mutual information rate, Ip, exchanged between 
the variables x and y is given by the rate of information 
produced in the one-dimensional space of the variable x, 
denoted as Hx, plus the rate of information produced in 
the one-dimensional space of the variable y, denoted as 
Hy , minus the rate of information production in the (x, y) 
space, denoted as H^y But, H^ = Hy = max(Ai,A2), 
and IIxy=IlKs=^i + A2, if (Ai, A2) > 0, Hxy—HKS—^i, 
if A2 < 0, and IIxy=IlKS=M^ otherwise. 

So, either Ip = Ai — A2, case that happens for when 
s — -fl, or Ip = A2 — Ai, case that happens for when 
s — —1. Ifs = -|-1, the larger the coupling strength, 
the smaller the KS-entropy, Hks- If s ~ ~1: the larger 
the coupling strength, the larger Hks- In fact, as we 
discuss further, Eq. P^ for s = —1 is a model for a 
self-excitable channel, and for s = -1-1 is a model for a 
non-self-excitable channel. In either case, the MIR can 
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be calculated by using the modulus operation as in Ip = 
|Ai — A2I. For larger networks, one can generalize such 
an equation using the conditional exponents arriving to 
an equation of the form as presented in Eq. ^. 

This equation points out to a surprising fact. Even 
when the level of desynchronization in Eq. (fT4|) is larger 
(A^ > A""^), which happens when s = — 1, there is a posi- 
tive amount of information being transferred between the 
two variables. 

In Figure [TOl we show the phase space of Eq. (fH)) 
for a coupling strength equal to (7=0.237. In (A), we 
illustrate a typical situation that happens in non-self- 
excitable channels (s = -1-1). The elements become 
synchronous presenting a trajectory that most of the 
time lies on the synchronization manifold defined by 
Xn — yn=0. In (B), we show a typical situation that hap- 
pens in self-excitable channels (s = —1). The elements 
become non-synchronous presenting a trajectory that lies 
on the transversal manifold defined by Xn+yn — c—0, with 
c being a constant within the interval c € \0. 1]. 
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FIG. 10: Trajectory of Eq. ^ for o-=0.237 with s 
(A) and s = 1 in (B). 
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In (A) , by observing the variable a;„ one can correctly 
guess the value of y„ since a;„ = y„. Apparently, that 
is not the case in (B): by observing the variable a;„, one 
might have difficulty in guessing the value of the variable 
y„, since c G [0,1]. Notice that the larger the amount 
of information being exchanged between Xn and ?/„, the 
larger the chance that we guess correctly. In order to 
estimate the amount of information being exchanged be- 
tween Xn and y„, we proceed in the following way. 

For the non-self-excitable channel (s=-|-l), we coarse- 
grain the phase space in I? small squares. Each square 
has one side that represents an interval of the domain 
of the variable a;„ and another side which is an in- 
terval of the domain of the variable 2/„- Calling p; 



the probability that a trajectory point visits the inter- 
val Xn = [(« — l)/i,i/-^]; with i — 1,...,L, and py , 
the probability that a trajectory point visits the inter- 
val y„ — \(i — 1)/L,i/L], and finally, Px-'y , the prob- 
ability that a trajectory point visits a square defined 
by Xn = [{i - 1)/L,i/L], y„ = [{j - l)/L,j/L], with 
j = 1, . . . , L, then, the MIR between x„ and j/„, denoted 
by /, is provided by 

/ = -l/log(L)[-^log(p('))-^log(pW)+^log(p(?;,^')]. 

i i i,j 

(15) 
Notice that the evaluation of the MIR by Eq. p^ under- 
estimates the real value for the MIR, since Eq. p^ is a 
dynamical system and the information produced by the 
dynamical variables (for example the term — ^ • log {px ) 
that measures the information produced by the variable 
Xn) should be provided by conditional probabilities, i.e., 
the probability that a trajectory point has of visiting a 
given interval followed by another interval, and so on, in 
fact the assumption used to derive Eq. ([5]). 

In Fig. [TlTA). we show the phase space of Eq. (fH|) 
with s=-|-l and for a = 0.237. In Fig. [TlTB). we show 
by the plus symbol, Ip, as calculated by Eq. ([5]) and by 
circles, /, as estimated by Eq. ([T5|. 
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FIG. 11: Results for Eq. CH) with (tO.237 and s=-l [shown 
in (A) and (B)] and for s=-l [shown in (C) and (D)]. Phase 
space of Eq. (|14|l and in (B), the MIR as calculated by Eq. 
([5]) and as estimated by Eq. (|15p . (C) Phase space of Eq. 
(|14p in the new coordinate frame XnVS.Xn+i and in (C), the 
MIR as calculated by Eq. ([5]) and as estimated by Eq. (|16p . 
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For the self-excitable channel (s=-l) Eq. p3|) sup- 
plies a null MIR, and therefore, it can no longer be used. 
But, as discussed in [ijj, the MIR can be coordinate 
dependent, and one desires to have the coordinate that 
maximizes the MIR. Aiming at maximizing the MIR, 
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when the channel is of the self-excitable type, we trans- 
form Eq. (|14p into an appropriate coordinate system, 
along the transversal manifold, where most of the in- 
formation about the trajectory position is located. We 
define the new coordinate as X„ = 1/2 (a;„ — j/n -I- 1) and 
X„+i = l/2(a;„+i-y„+i + l). The trajectory (X„,X„+i) 
in this new coordinate system [for the same parameters 
as in Fig. MB)] is depicted in Fig. [IHC). 

The MIR being transferred between X„ and Xn+i is 
related to the knowledge we acquire about X^+i by ob- 
serving Xn, or vice- versa. In Fig. [TlTC). we can only 
be certain about the value of Xn+i, when Xn is close to 
either or 1. 

To estimate the MIR, we recall that an encoded version 
of such a dynamical system can be treated as a symmetric 
binary channel of communication. Xn is regarded as the 
transmiter and X„+i is regarded as the receiver. When- 
ever the map in the transformed coordinates XnVS.Xn+i 
is non-invertible, we consider that by making measures of 
the trajectory point Xn+i one cannot guarantee the ex- 
act position of the trajectory of X„, which constitutes an 
error in the transmission of information. Whenever the 
map is invertible, by measuring the trajectory of X„-|_i 
one can surely know the exact position of the trajectory 
Xn , which corresponds to a correct transmission of infor- 
mation. Calling, p the probability at which the map is 
invertible, then, the MIR between Xn and Xn+i is given 
by 



/e = 1 + (1 - P) log {l-p)+p log (p). 



(16) 



The value of lO/g for Eq. pA^ with s = — 1 are shown 
in Fig. [TIT D) by circles. The theoretical value, Ip, pro- 
vided by Eq. ([5]) is shown by the plus symbol. 

A final comment on the characteristics of a self- 
excitable channel and of a non-self-excitable channel is 
that while in a self-excitable channel the larger the syn- 
chronization level, the larger the MIR but the smaller 
the KS-entropy, in a non-self-excitable channel the larger 
the desynchronization level, the larger the MIR and 
the larger the KS-entropy. Note that i//i-5=21og(2) + 
log(l-2scr), for cr < 0.25. 



B. Positiveness of the MIR for self-excitable 
channels in the (non-linear) HR network 

To show that indeed Ip should be positive in case of a 
self-excitable channel in the HR network, one can imagine 
that in Eq. ^ the coupling strength is arbitrarily small 
and that N~2. At this situation, the Lyapunov exponent 
spectra obtained from Eq. ([3]) are a first-order perturba- 
tive version of the conditional exponents, and they ap- 
pear organized by their strengths. One arrives at Ai = A^ 
and A2 = A^, which means that the largest Lyapunov ex- 
ponent equals the transversal conditional exponent and 
the second largest Lyapunov exponent equals the condi- 
tional exponent associated with the synchronous mani- 
fold, i.e., the Lyapunov exponent of Eq. Q. Using sim- 



ilar arguments to the ones in Refs. [llj, OJ, |29| , we have 
that the MIR is given by the largest Lyapunov exponent 
minus the second largest, and therefore, /c = Ai — A2, 
which can be put in terms of conditional exponents as 
Jp <A2-Ai. 



C. The inequality in Eq. (O 

To explain the reason of the inequality in Eq. ([5]), 
consider the nonlinear term in Eq. p4p is non null and 
s—l, and proceeds as further. 

For two coupled systems, the MIR can be writen in 
terms of Lyapunov Exponents [13, [30| . For two coupled 
systems, the MIR can be exactly calculated by Ic = 
Ai — A2, since A" = Ai and A^ = A2, assuming that 
both Ai and A2 are positive. Calculating the conditional 
exponents numerically, we can show that Ip > Ici E^nd 
thus Ip is an upper bound for the MIR. For more details 
on this inequality, see [3ll | 



XVII. BUST PHASE SYNCHRONIZATION 

(BPS) 

Phase synchronization [32| is a phenomenon defined by 



\A(l,{kJ)\ = \(j>k - m^i\ < r, 



(17) 



where cjjk and 0/ are the phases of two elements Sk and 
Si, m = ioi/ujk is a real number [33], where lu^ and uji are 
the average frequencies of oscillation of the elements Sk 
and Si, and r is a finite, real number [3^. In this work, 
we have used in Eq. (J17p m = 1, which means that we 
search for oj^ ■ uJi~l:l (rational) phase synchronization 
32]. If another type oiujk '■ ijJi-PS is present, the methods 
in Refs. [3J, [Sa, ISa] can detect it. 

The phase 4> is & function constructed on a 2D sub- 
space, whose trajectory projection has proper rotation, 
i.e, it rotates around a well defined center of rotation. So, 
the phase is a function of a subspace. Usually, a good 2D 
subspace of the HR neurons is formed by the variables 
X and y, and whenever there is proper rotation in this 
subspace a phase can be calculated as shown in Ref. [37| 
by 



M 



yx — xy 



y 



dt. 



ix" 



(18) 



If there is no proper rotation in the subspace (x, y) one 
can still find proper rotation in the velocity subspace 
(i, y) and a corresponding phase that measures the dis- 
placement of the tangent vector [33 can be calculated as 
shown in Ref. [33| by 



^,{t) 



yx — xy 

(±2+2/^) 



dt. 



(19) 



If a good 2D subspace can be found, one can also define 
a phase by means of Hilbert transform, which basically 
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transforms an oscillatory scalar signal into a two compo- 
nent signal [3^]. In the active network of Eqs. ([7]) with 
an all-to-all topology and N=4:, for the coupling strength 
interval a = [0,0.05], the subspace (x,y) has proper ro- 
tation, and therefore, (f>s{t) is well defined and can be 
calculated by Eq. (fT8|) . However, for this coupling inter- 
val, Eq. ([T7]) is not satisfied, and therefore, there is no 
PS between any pair of neurons in the subspace {x,y). 

For the coupling strength interval a = [0.05, 0.24], the 
neurons trajectories lose proper rotation both in the sub- 
spaces {x,y) and {x,y). In such a case, neither 0i(i) nor 
4>v{t) can be calculated. This is due to the fact that the 
chaotic trajectory gets arbitrarily close to the neighbor- 
hood of the equilibrium point (x,y)=(0,0), a manifesta- 
tion that a homoclinic orbit to this point exists. 

In fact, the Hilbert transform fails to provide the phase 
from either scalar signals x or y, since these signals do 
not present any longer an oscillatory behavior close to 
the equilibrium point. In such cases, even the traditional 
technique to detect PS by defining the phase as a func- 
tion that grows by 27r, whenever a trajectory component 
crosses a threshold cannot be used. Since the trajec- 
tory comes arbitrarily close to the equilibrium point, no 
threshold can be defined such that the phase difference 
between pairs of neurons is bounded. Notice that by 
this definition the phase difference equals 27rA7V, where 
AN is the difference between the number of times the 
trajectory of Sk and Si cross the threshold. For the neu- 
ral networks, AA^ could represent the difference between 
the number of spikes between two neurons. A spike is 
assumed to happen in Sk if Xk becomes positive. 

In order to check if indeed PS exists in at least one 
subspace, alternative methods of detection must be em- 
ployed as proposed in Refs. [33, [Sg]- In short, if PS 
exists in a subspace then by observing one neuron tra- 
jectory at the time the other bursts or spikes (or any 
typical event), there exists at least one special curve, P, 
in this subspace, for which the points obtained from these 
conditional observations do not visit its neighborhood. A 
curve P is defined in the following way. Given a point xq 
in the attractor projected onto the subspace of one neu- 
ron where the phase is defined, P is the union of all points 
for which the phase, calculated from this initial point xq 
reaches n{r), with n = l,2,3,...,c» and (r) a constant, 
usually 2tt. Clearly an infinite number of curves P can be 
defined. For coupled systems with sufficiently close pa- 
rameters that present in some subspace proper rotation, 
if the points obtained from the conditional observations 
do not visit the whole attractor projection on this sub- 
space, one can always find a curve P that is far away from 
the conditional observations. Therefore, for such cases, 
to state the existence of PS one just has to check if the 
conditional observations are localized with respect to the 
attractor projection on the subspace where the phase is 
calculated. 

Conditional observations of the neuron trajectory Sk in 
the subspace {x,y), whenever another neuron Si spikes, 
in the system modeled by Eqs. ([7]) with a star coupling 
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FIG. 12: The network of Eqs. (0 with a star configuration with 
N=A, and a=0.265. The curve F, a continuous curve transversal to 
the trajectory, is pictorially represented by the straight Hne F. (A) 
the green line represents the attractor projection on the subspace 
{x, y) of the neuron S2, and red circles represent the points obtained 
from the conditional observations of the neuron 52 whenever the 
neuron 54 spikes. The point (x, y) = (0.0) does not belong to F. 
(B) Green dots represent the reconstructed attractor 22(4) X Z2{t — 
t), for T=30, and red circles represent the points obtained from the 
conditional observation of neuron 52 , whenever the reconstructed 
trajectory of the neuron 54 crosses the threshold line 24(4 — r) = 
3.25 and Z4(t) > 3. 



topology and N=A, are not localized with respect to a 
curve P, for the coupling strength a < dps- An exam- 
ple can be seen in Fig. [TWA), for a = 0.265. The set 
of points produced by the conditional observations are 
represented by red circles, and the attractor by the green 
points. Therefore, there is no PS in the subspace {x,y). 

In order to know on which subspace PS occurs, we 
proceed in the following way. We reconstruct the neuron 
attractors by means of the time-delay technique, using 
the variable z. This variable describes the slow time- 
scale, responsible for the occurrence of bursts. The re- 
constructed attractor z{t) x z(t — r) has proper rotation 
[see Fig. [T^B)] and the points obtained from the con- 
ditional observations do not visit the neighborhood of a 
curve P, then, there is PS in this subspace. Indeed, we 
find localized sets with respect to a curve P in the recon- 
structed subspace (z(t) x z{t — t)), for a > 0.265. So, 
ctbp5=0.265. 

So, for the coupling a — [aBPS,o-ps[i there is no PS 
in the subspace {x,y) but there is PS in the subspace of 
the variable z. In this type of synchronous behavior, the 
bursts are phase synchronized while the spikes are not. 
This behavior is regarded as bursting phase synchroniza- 
tion (BPS). For simplicity in the analyses, we say that 
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BPS happens when for at least one pair of neurons there 
is phase synchronization in the bursts. Phase synchro- 
nization (PS) happens in the network when the average 
absolute phase difference 



-^^|A0,(fc,O|, 



with k — l,N —1 and / ~ k + 1, N among all the pairs of 
elements, is smaller than 2tt, with the phases defined by 
either Eq. (fT8|) or Eq. (fT9|) , where the index L represents 
either the index s or v. Furt her, we say complete syn- 
chronization (CS) takes place [i3|, when the variables of 
one neuron equal the variables of all the other neurons. 

For the analyses in this work, (Jbps represents the cou- 
pling parameter for which BPS first appears, i.e., BPS 
exists if cr > ubps ■ <^ps represents the coupling parame- 
ter for which PS first appears, i.e., PS exists if cr > ups- 
Finally, acs represents the coupling parameter for which 
CS first appears, i.e., CS exists if cr > cfcs- There might 
exist particular parameters for which PS (or BPS) is lost 
even if cr > aps (resp. a > ubps)- But these parameters 
are not typical and we will ignore them. For example, in 
the network composed by 6 elements with the nearest- 
neighbor topology [Fig. EIJB)], for a = 0.825 PS is lost. 

Note that these phenomena happen in a hierarchical 
way organized by the " intensity" of synchronization. The 
presence of a stronger type of synchronization implies in 
the presence of other softer types of synchronization in 
the following order: CS ^ PS ^ BPS. 



and delete all links connected to it. In the following, we 
choose randomly a new degree k to this element and con- 
nect this element (in a bidirectional way) to k other ele- 
ments randomly chosen. This procedure generates a new 
network that possesses the Laplacian Q' , whose eigenval- 
ues produce a value B'. To decide if this mutation is 
accepted or not, we calculate Ae = B' — Bo- If Ae > 0, 
the new network whose Laplacian is G' is accepted. If, 
on the other hand, Ae < 0, we still accept the new mu- 
tation, but with a probability p(Ae) = exp(— A/eT). If a 
mutation is accepted then the network whose Laplacian 
is Qo is replaced by the network whose Laplacian is Q' . 

The parameter T is a kind of "temperature" which con- 
trols the level of noise responsible for the mutations. It 
controls whether the evolution process converges or not. 
Usually, for high temperatures one expects the evolution 
never to converge, since new mutations that maximizes 
B are often not accepted. In our simulations, we have 
used T ^ 0.0005. 

These steps are applied iteratively up to the point 
when |Ae| = for about 10,000 steps, being that we con- 
sider an evolution time of the order of 1,000,000 steps. 
That means that the evolution process has converged af- 
ter the elapse of some time to an equilibrium state. If 
for more than one network topology |Ae| = for about 
10,000 steps, we choose the network that has the larger 
B value. 

This constraint avoids the task of finding the most op- 
timal network topology. However, we consider that a 
reasonably low number of mutations would recreate what 
usually happens in real networks. 



XVIII. EVOLUTIONARY CONSTRUCTION OF 
A NETWORK 

In our simulations, we have evolved networks of equal 
bidirectional couplings [331 ■ That means that the Lapla- 
cian in Eq. ([1]) is a symmetric matrix of dimension N 
with integer entries {0, 1} for the off diagonal elements, 
and the diagonal elements equal to — X]i Sij, with i ^ j. 

Finding the network topologies which maximize B in 
Eq. ([5]) is impractical even for moderately large N. Fig- 
uring out by " brute force" which Laplacian produces the 
desired eigenvalue spectra would require the inspection of 

a number of j^j. configurations. To overcome this 

difficulty, Ref. 14] proposed an evolutionary procedure 
in order to reconstruct the network in order to maximize 
some cost function. Their procedure has two main steps 
regarded as mutation and selection. The mutation steps 
correspond to a random modification of the pattern of 
connections. The selection steps consist in accepting or 
rejecting the mutated network, in accordance with the 
criterion of maximization of the cost function B, in Eq. 

We consider a random initial network configuration, 
with A'^ elements, which produce an initial Laplacian Qq, 
whose eigenvalues produce a value Bq for the cost func- 
tion. We take at random one element of this network 



XIX. CONSTRUCTING A NETWORK FROM A 
SET OF EIGENVALUES 

Given a N x N Laplacian matrix G, we can diagonalize 
it by an orthogonal transformation, viz 



0^.g.0 = 7l, 



(20) 



where 1 represents the Unity matrix, 7 represents the 
vector that contains the set of eigenvalues 7; oi Q {i = 
l,...,N), and O is an orthogonal matrix, O.O-^ = 
O-'^.O = 1, whose columns are constructed with the or- 
thogonal eigenvectors of Q, namely O = [x/i, ■i/2, . . . , utv]- 
Accordingly, 



g^o.ji.o^, 



(21) 



which means that Q can be decomposed into a multipli- 
cation of orthogonal matrices. By using the spectral form 
of Eq. (I^T]) , the Laplacian Q can be calculated from 



N 



Q = ^Vi.-li.vJ . 



i=l 



(22) 



Any other Laplacian, Q' ^ can be constructed by using 



19 



the set of eigenvalues 7, viz 



N 



<^. 



(23) 



Of course, in order for the active network that is con- 
structed using Q' to present the synchronization man- 
ifold Xi = X2 — x^ = ■■■,— a;„, the vector v[, with 
N elements, is given by v'^ — -^[1, 1, 1, 1, ■ • ■ , 1], and 
the other vectors are found by choosing arbitrary vectors 
v[ which are made orthogonal using the Gram-Schmidt 
technique. 
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